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The van der [Nil limit eveles arc generated at small amplitudes hv the computer implementation 
of the Poineare-Lindstedl method. The formal algebraic solution is accomplished by manipula- 
tions ol Poisson scries, and the FORTRAN programming of the inductive algorithm yields the 
phase-shifting limit cycles to graphical accuracy over the range "S A. < 1.5. This improves 
upon the method of Deprit and Horn in two ways. First, because the formal solution is carried 
out by hand, an algebraic processor is not necessary. Second, the standard solutions which 
they generated arc only valid for < A. < 1.2 whereas the phase-shifting limit cycles are still 
valid at A = 1.5; that is. they do not exhibit the Gibbs phenomenon even at A = 1.5. 
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1 . Introduction 

The computer generation of the Poincare-Lindstedt solutions to nonlinear differential equations is reported 
by Deprit and Rom in [1], and the term Poisson series for combined power-Fourier series is coined in |2|. 
More recently the Poincare-Lindstedt method has been reapplied to different forms of the nonlinear oscillator 
equation in [3], [4], and [5]. 

In this paper the lead of [6] is followed, and the limit cycles of the van der Pol equation are constructed at 
small amplitudes. In the method of Deprit and Rom, a machine language algebraic processor is used to 
manipulate the Poisson series. It is shown here that it is not necessary to use such involved software tools 
because these algebraic operations can be carried out by hand, and the resulting system of recursion relations 
can be implemented in a higher level language such as FORTRAN. The bonus for first explicitly developing 
the recurrence relations for the van der Pol equation is that the existence of the phase-shifting cycles becomes 
apparent, and they can be generated to a larger value of the parameter than the standard cycles. 



The van der Pol equation 



2. Solution by Poisson Series 



Y"(t) - \(1 -Y 2 )Y' +Y =0, (1) 



has enjoyed considerable interest in the literature of applied mathematics because it has singular, periodic 
solutions called limit cycles. The approach of the general solution to the limit cycle can be seen in the results 
of Davis in [7]. The cycles have been computed for all values of the parameter by purely numerical methods 
by Urabe in [8) and Clenshaw in [9], In fact the numerical production of cycles is so simple it is given as an 
excercise in [ 1()|. 

A more difficult task is their analytical production. The method of Poincare-Lindstedt as applied to (1) is 
well known and can be found in Cesari [11]. As noted in [2] and [3], an advantage to the computer 
implementation of the Poincare-Lindstedt method is the speed and accuracy with which the whole process can 
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be performed. That is, within the interval of convergence the numerical evaluation of the Poisson series by 
Horner's rule and Goertzel's algorithm, see Coffey and Deprit [14], yields numerical values very much more 
rapidly than through normal methods of numerical integration. 

The method of [5] is extended to the van der Pol equation by the assumption that the limit cycles have the 
Poisson series expansion (1.3) in table 1, where i = (— 1) 1/2 . The frequency of the cycle co is defined by r 



TABLE 1. Definitions of Coefficients 



(t) — 



2 X 2 '"" 2 co m , (1.1); o> 2 = 2 k 2m -m m , (1.2) 

m = \ m = \ 

00 00 

¥--■■ 2 ««*-»" 2 k* m -*(R\»j + i\I\»j), (1.3) 

j=-oo m=l 

00 00 

Y 2 = £ ,OHIr 2 \2m-2 (/? «2). + fc/g^), (1.4) 

j=-oo m=\ 

\Y 2 Y' = 2 e l(2J - 1)T 2 ^ 2m_2 (RS!j + iX/Sjj), (1.5) 

j=— ao m=l 

= co t where the r-period is 2 n. The algorithm is simplified by expanding both co and to 2 as given in (1.1) 
and (1.2). By standard series manipulations, the problem of integrating (1) is reduced to the inductive 
algorithm of table 2. The indicial conditions on the algorithm are 

R[\\ = *>, =n, = 1. (2) 

The reflection relations are also given in table 2; they are used to obtain real solutions and to limit storage of 
coefficients to positive values of/. At the end of each equation in table 2 the selection rules are given which 
each nonzero coefficient must satisfy for positive/. For indices outside of the range of the selection rules, the 
corresponding coefficient is zero. 

In table 2 recursion relations (2.1) and (2.5) are derived by substitution of (1.3) in the left-hand side of 
(1.4) in table 1 and by equating coefficients of equal powers of \e 17 . Recursion relations (2.2) and (2.6) are 
obtained by substitution of (1.1), (1.3), and (1.4) in the left-hand side of (1.5). Because the nonlinear term in 
(1) can be expressed algebraically as the Poisson series (1.5), which is of the same form as the solution series 
(1.3), it might be called the nonlinear closure property of Poisson series. 

By substitution of (1.3) and (1.5) in (1), two solution relations are obtained as 

m rn-l 

- (2; - 1 ) 2 1 ft m _ k+i R gj + (2; - 1) 2 « m-kl B + R Z + R l A = 0, (3 ) 

and 

m m 

- (2; - i ) 2 2 n ,„-, +1 / B + (2/ - 1) 2 w »-* +1 R K + 1 Z + I SJ'j = o. (4) 

k=\ k=\ 

When/ > 2, relation (3) is solved with the indicial conditions (2) as (2.9). When / = 1, relation (3) gives by 
(2.4) the coefficients in the square of the frequency as (2.7). This is Lindstedt's method of casting out secular 
terms. Fory = 1, relation (3) leaves arbitrary the sequence of linear amplitudes R^ . 

Wheny > 2, relation (4) is solved by (2) as (2.3). When j — 1, relation (4) leaves I^li arbitrary, but by 
(2.4) the sequence of linear amplitudes is determined by the algebraic equation 
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Table 2. Algorithm for Phase-Shifting Limit Cycles 



i<«.-i»>.£ c2,-i>r£»< 2 > , . ,E-i ,* (1) -fi <2 L- ,!>» ' (1) 1 

m,i m,l-j *-• K /-. m-k+1 , j-p+1 *-*. k-n+1 n,p /-< m-k+2 , i-p+1 ^, k-n n,p 

J J p=-«o L k=l J n=l k=2 J r n*l -I 

m,j m,l-j |_ m,j ^ m-k+1 k,j £^ m-k+1 k,jj 



The m-th approximation is complete, replace m by m + 1. 



!«>..,<« . . £ £(,<» . /".I«> . ,5") 

m,j m,2-j y?m£m\ m-k 'J"P k 'P m-k,j-pTi f p 



R m,j " Vl-J = "2 <*"*> £* k+1(J _ p+1 L " k _ n I n , P + ^m-k+l.j-p+l^k-n + l^.p » 
J J p=-e» L k=2 J r n*l r k=l J r n=l J 



m m, 1 , z — ' m-k + 1 k, 1 

k=2 



m,j m,l-j [ m,j ££ m-k k,j J £* " m-k+l~k, J J 



Set R . ■ and compute (2.1) and 12.2), then compute 



i cl) - k\*T , r (1) i (3) 1 
R m,i *[^» l w «-k*l , k f l «,1 J ' 



Return to (2.1). 



2m > j > 1 . (2.1) 



2m > j > 1 . (2.2) 



(j-1) , 2m > j > 2. (2.3) 



l il ] - 0. (2.4) 

m.l 



2m-l > j > 1. (2.5) 



2m-l > j > 1. (2.6) 



(2.8) 



/4j( j-1) , 2m- 1 > j > 2. (2.9) 



/ ( ™ 3 !i-2 oj m _ k+ M\ = o. (5) 

fc=l 

The coefficients R j^j occur explicitly ill the second term in (5) and implicitly in the first. To make this 
dependence explicit, (2.1) is rewritten as 

RZ = 4 Rj&, + R{&, mdRMj = 2 IV,!,', + «<&, (6) 

where the remainder RjIEj is independent of R^Ji. Relation (2.2) is rewritten by (6) as 

/',;::, = 3 «& + /<»„ (7) 
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where the remainder 1 { ^ A is also independent of Rmli* The remainder terms are easily found in a computer 
program by setting Rm\ to zero and by then executing the subroutines for (2.1) and (2.2). Substitution of (7) 
in (5) gives (2.10) upon solution for R { ,n A . Thus, the dissipation terms in the van der Pol equation determine 
the sequence of linear amplitudes. 

3. Standard and Phase-Shifting Cycles 

When j = 1, relation (4) leaves arbitrary the phase amplitudes /{JJj and unlike R ( m,i their values are not 
determined by the van der Pol equation itself. This harks back to the fact that (1) is autonomous. In [6], [7], 
[8], and [9] the phase of the cycle is determined by imposing the initial condition 

F'(0)=0, (8) 

for all values of A. That is, the phase of the cycles is adjusted so that the initial amplitude is a maximum. To 
fulfill this condition, it is found from (1.3) that the phase amplitudes must satisfy 

2?fl 

/£. = "2 (2/-D/&. (9) 

3=2 

Because so many investigators impose the condition (8), here the resulting solutions are called the standard 
cycles. 

The arbitrary sequence Im]i can De chosen, however, in other ways. In particular the simple choice (2.4) 
can be made. When /£J*j is not chosen to satisfy (9), then the condition (8) no longer holds. The resulting 
cycles are said to be phase-shifting, and table 2 is arranged to generate coefficients for phase-shifting cycles. 
The algorithm for the standard cycles is obtained by replacement of (2.4) by (9), by replacement of (2.7) by 

m—\ m—\ 

\l m = R\n,i ~~ 2* *Lm-k+\Rk,l + 2* <*> m -klk,n 
k=2 fc=l 



and bv replacement of (2.10) bv 



(tom-k+iRkA + **m-k+Jk,l) J 



(3) 



These replacements require the algorithm for the standard solutions to have more arithmetic operations than 
for the phase-shifting solutions. It is probably for this reason that the phase-shifting cycles can be produced 
to larger values of the parameter than the standard cycles. 

4. Comparison of Results 

To compare results with previous investigators, the standard cycles were generated. The amplitude of the 
standard solutions is 

00 

a(k) -r(X,0) = X a m k 2 '"-\ (10) 

m = \ 

where the amplitude coefficients are computed from the linear amplitudes by 

2m- 1 
a = 2 y R (1) 

Other than by direct comparison, there is no simple way to determine how the various sources of error affect 
the final results. To study the propagation of error, the program was executed on an IBM 360/75 at the 
University of Illinois at Urbana-Champaign in single precision which gives six to seven decimal digits. The 
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first column of table 3 contains the results of Deprit and Rom who used double precision on an IBM 7094, 
and the second column contains results of this investigation. The coefficients in the two columns show perfect 



agreement with each other and with the rational fracti 



l)\ fiand by Clenshaw in [9] only to order m 



i\\\ two decimal agreement. That 
man\ significant digits. What is 



= 3. At order m — 4 there is a disadvantageous subtraction, and then 

there is such a disagreement is expected since Deprit and Rom used 

surprising is that the disagreement occurs at such a low order. The fourth order is beyond my ability for hand 

computation, and it would be interesting to sec the method of |4| applied to (I) to give at least the m = 4 

coefficients as rational fractions. 



TABLE 3. Amplitude and Frequency Coefficients 



m 


a m ( Deprit & Rom) 


a m 


<o m 


r TO 


1 


2 


2 


1 


1 


2 


0.010 416 666 667 


1.04167E-2 


-6.25000E-2 


2.52 


3 


-0.00] 868 127 893 


-1.86813E-3 


5.53385E-3 


2.83 


4 


0.000 018 204 203 


L.83739E-5 


3.95598E-5 


4.26 


5 


0.000 062 r>2() 749 


6.25358E-6 


L.33220E-4 


2.70 


(> 


o.ooo 007 o8o 633 


-7.98043E-6 


L.22798E-5 


2.80 


7 


o.ooo 002 189 064 


-2.22516E-6 


4.52331E-6 


2.58 


8 


o.ooo ooo 678 802 


6.76593E-7 


L.09645E-6 


2.50 


9 


0.000 000 048 985 


4.81233E-8 


-1.25366E-7 


2.55 


10 


-0.000 000 046 470 


-4.29809E-8 


7.70256E-8 


2.37 


II 


o.ooo 000 002 256 


2.20042E-9 


-1.04284E-9 


2.68 


12 


0.000 000 002 732 


2.90587E-9 


-4.70514E-9 


2.30 


13 


-0.000 000 000 465 


-8.84975E-10 


6.03579E-10 


2.34 


14 


-0.000 000 000 130 


-1.49888E-10 


2.41538E-10 


2.27 


L5 


0.000 000 000 047 


3.97400E-] 1 


-6.86656E-] 1 


2.24 


L6 


0.000 000 000 042 


5.08794E-] 1 


-8.57589E-12 


2.28 


17 




-5.45275E-12 


5.78845E-12 


2.10 


18 






-7.64649E-14 


2.37 


19 






-4.04249E-13 


2.16 


20 






5.35501 K- 14 


2.19 


21 






2.30394E-14 


2.15 


22 






-6.68763E-15 


2.14 



The agreement improves slightly at /// = 5 and 6, but it worsens as higher orders are reached. On the other 
hand, the Poincare-Lindstedt method is numerically stable because even at order 15 the amplitude' 
coefficients disagree with a relative error of only 15 percent. Although Deprit and Rom report numerical 
values to order 16, there is an additional unknown correction in their last coefficients. Evidently they output 
their coefficients after (2.9) but before evaluating (2.10). A better output point is after (2.4) as then the order 
m is completed in (11) and (12). 

The disagreement in the amplitude coefficients is not significant, however, because it is the values of a{k), 
o>(X), K(X, r) and Y'(t) which are more important. The values of the amplitude and frequency are given in 
table 4. The purely numerical results of [6] or [9] are listed in the first and third columns. The IBM single 
precision evaluation in the second column shows perfect agreement for < X < 1.5 and is in relative error of 
— 0.02 percent at X = 1.75. The value obtained by Deprit and Rom for a (1.75) is slightly better with a 
relative error of —0.004 percent. 

Even though Deprit and Rom used twice as main digits in their implementation of the Poincare-Lindstedt 
method, they were unable to produce limit cycles without the Gibbs phenomenon beyond k = 1.2. The reason 
for this is apparent from figures 1.1 and 1.2. The imposition of the standard condition (8) apparently causes 
the early onset of the uncontrolled oscillation. In figure 1. 1 the standard cycles are plotted for order 17 for X 
= 0, 0.5, 1.0, and 1.5 as obtained in IBM single precision. 

Figures 1.1, 1.2, and L.3 overlap so that figure 1.1 refers to the first two families of curves. Figure 1.2 
refers to the second pair of families. The top family of curves are F(X, r) for X = 0, 0.5, 1.0., and 1.5; and the 

597 




FIGURE 1. 1.1 shows standard cycles; 1.2 shows phase- 
shifting cycles; 1 .3 shows roots. 



second family of curves are Y'(t) for the same values of X. The curves of figure 1. 1 represent the standard 
cycles computed by the method of this paper. The dashed lines for Y'(t) in figure 1.1 are obtained from the 
m = 16- and 17-th order polynomials and clearly exhibit the Gibbs phenomenon at X = 1.5. The onset of the 
Gibbs phenomenon occurred in Deprit and Rom's results at X = 1.35. 

In figure 1.2, the first family of curves are the phase-shifting limit cycle solutions K(X, r) for X = 0, 0.5, 
1.0, and 1.5. The second family in figure 1.2 is Y'(t) for X = 0, 0.5, 1.0, 1.5, and 1.75. Thus, it is possible 
to generate the phase-shifting limit cycles of the van der Pol equation without the Gibbs phenomenon at X = 
1.5. The Gibbs phenomenon is, however, exhibited at X = 1.75 in the phase-shifting cycles as can be seen 
in figure 1.2 where the m = 21-st and 22-nd order polynomials are plotted as dashed curves. 

5. The Phase-Shifting Cycles 

When the algorithm of table 2 was executed in IBM single precision, it was found that unlike the standard 
solutions there is a negligible amount of the Gibbs phenomenon in the phase-shifting solutions at X = 1.5 for 
m = 15. The program for the phase-shifting cycles was again run on a CYBER 175 where the program 
achieved the 22nd order in 130 seconds. The 60 bit, single precision word of the CYBER accumulator gives 
14 decimal digits of significance. All of the frequency coefficients are listed to six decimals in the third 
column of table 3, and the solution coefficients are given to the 6th order in tables 5 and 6. In the fourth 
column of table 4 values of the frequency are given, and there is perfect agreement with the third column for 
< X < 1.5. The 22nd order phase-shifting cycles are plotted in figure 1.2. In addition, Y '(t) is plotted for 
X = 1.75 to show the Gibbs phenomenon at orders 21 and 22 as dashed lines. Thus, by not imposing the 
constraint (8), the Poincare-Lindstedt method of constructing the limit cycles of the van der Pol equation 
yields at least graphically accurate results over the interval < X < 1.5. 

6. Convergence of the Series 

Although the Poisson series solution to (1) is shown rigorously to converge at very small values of X by Hale 
in [12], it is not obvious for how large a value of the parameter the series continues to converge. Because the 
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Tablk 4. Amplitude and Frequency 






k 


a(\), [6, 9] 


a(X), m = 17 


w(\), [6, 9] 


o>{ A). /// = 22 


(t>{K). m = 3 


0.5 


2.002 487 9 


2.002 487 


0.984 721 


0.984 721 


0.984 721 


1.0 


2.008 619 ( > 


2.008 619 


0.942 955 8 


0.942 956 


0.943 034 


1.5 


2.015 22(> 5 


2.015 227 


0.885 407 ( > 


0.885 408 


0.887 390 


1.75 


2.017 860 1 


2.018 208 


0.854 450 5 


0.854 430 


0.860 495 



character of the solution to the linearized form of (1) changes at X = 2 from underdamped to overdamped, it 
is expected that X = 2 is the radius of convergence of the series in table 1. When scries solutions arc 
generated on a computer, the radius of convergence is not usually known beforehand; and it is helpful if the 
program can determine the range of validity of the series. The root test used in [13] is more reliable than the 
ratio test used in [5]. It can be seen in the last column of table 3 and in figure 1.3 that the roots 



appear to have the lower bound of 2. 

By comparing the accuracy of the approximating polynomials at order/// with the behavior of the roots, it is 
found that the most accurate polynomials are those for the order just before a large increase in the root. For 
example, in table 3 there is a jump in the root at /// = 4; and it is seen in the last column of table I that a good 
approximation to the frequency is obtained at the 3rd order. 



7. Lower Order Approximations 

The real form of the phase-shifting cycles is found from (1.3), the reflection and selection rules as 



Y(K t) = 2 X X 2 



2 RS& cos (2/ - 1 ) T - X 2 tmj sin (2; - I ) r 



3=2 



(ID 



Bv the chain ride the derivative is found to be 



Y'(t) =-2wJ X 2 ' 



2 (2/ - 1 )R\ f )] j sin (2/ - 1)t + xj? (2/ - 1)/«>j cos (2/ 



It 



i=2 



(12) 



The greatest problem in the numerical evaluation of the limit cycles is the Gibbs phenomenon in the derivative 
(12). To see what order is required to obtain graphical accuracy for a given X, a TUTOR language evaluation 
of (1.1), (11), and (12) was used on the PLATO IV system at the University of Illinois. In figure 2 the phase- 
shifting cycles are plotted for X = 0, 0.25, 0.5, 0.75, 1.0, 1.25, and 1.5 which are obtained at orders L, 3, 
3, 3, 6, 9, and 15. The PLATO storage is rather limited, but it could accommodate the 510 solution 
coefficients necessary for order 15. 

Because the total storage of table 2 is only 7750 coefficients to order 25, it is perhaps not necessary to use 
a large computer to evaluate the algorithm. Since table 2 was coded in FORTRAN, it is probably possible to 
code it in the BASIC language of minicomputers. Although storage requirements arc not prohibitive, computer 
time is a limitation. The execution time for each order appeals to grow exponentially. For example, it took 3 
minutes of IBM 360/75 time to reach order 15. but an additional two orders to /// = 17 required two minutes 
more. 

In tables 5 and (>, the real and imaginary solution coefficients are given to six figures to the 6th order. 
There are sufficient coefficients in these tables and table 3 to construct the solutions (1.1), (11), and (12) on 
a programmable calculator to graphical accuracy over the interval < X < 1.0. 

This investigation of the van der Pol equation was originally undertaken to gain understanding of the limit 
cycle phenomena. The problem of the limit cycle behavior of periodic variable stars is reported in | 15|. 
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Figure 2. Plato display; van der Pol Limit Cycles. 



Table 5. Real, Phase-Shifting Solution Coefficients R { „ 



m/j 


L, 6, 11 


2. 7 


3, 8 


4, 9 


5, 10 


1 
2 


1.0 
7.81250E-3 


-4.68750E-2 


-2.60417E-2 






3 


-2.33968E-4 


4.10970E-3 


8.43189E-3 


6.23463E-3 


1.48926E-3 


4 


-1.51937E-4 

-6.17692E-4 


4.24933E-4 
-9.6241 7E-5 


-6.45271E-4 


-1.71476E-3 


-1.49991E-3 


5 


2.42362E-5 


-1.71589E-4 


-1.90775E-4 


5.73677E-5 


3.38734E-4 




3.62885E-4 


1.97997E-4 


5.64688E-5 


6.55160E-6 




6 


3.98953E-6 


9.46432E-7 


4.82361E-5 


6.52477E-5 


1.05810E-5 




-6.21963E-5 


-8.53587E-5 


-5.77692E-5 


-2.27698E-5 


-4.96343E-6 




-4.59207E-7 











Table 6. Imaginary, Phase-Shifting Solution Coefficients PJ^j 



m/j 


2. 7. 12 


3,8 


4, 9 


5, 10 


6, 11 


1 


0.125 










2 


-1.46484E-2 


-1.84462E-2 


-6.07639E-3 






3 


3.30183E-4 


3.04987E-3 


3.76636E-3 


1.99198E-3 


3.75231E-4 


4 


3.33934E-4 
-1.87977E-4 


1.02407E-4 
-2.49947E-5 


-5.44143E-4 


-8.28345E-4 


-5.57686E-4 


5 


-4.74168E-5 


-1.22377E-4 


-6.87602E-5 


8.32629E-5 


1.79145E-4 




1.48482E-4 


6.79942E-5 


1.68027E-5 


1.72960E-6 




6 


-9.60169E-6 


9. 72957E-6 


3.40962E-5 


2.76165E-5 


-7.88881E-6 




-3.68516E-5 


-3.77886E-5 


-2.16231E-5 


-7.47641E-6 


-1.45772E-6 




-1.22484E-7 
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